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Abstract We present two new methods for linear elasticity that simultaneously yield stress and 
displacement approximations of optimal accuracy in both the mesh size h and polynomial degree p. 
This is achieved within the recently developed discontinuous Petrov-Galerkin (DPG) framework. 
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element interfaces. We study locking-free convergence properties and the interrelationships between 
the two DPG methods. 
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1 Introduction 

In this paper we propose a new method for numerically solving the system of equations describing 
linear elasticity. The accurate computation of stresses is of critical importance in many applica- 
tions. Yet, many traditional methods only yield approximations to the displacement. This means 
that stress approximations must be recovered afterward by numerical differentiation. There are 
newer methods, of the mixed and discontinuous Galerkin (DG) category, which do give direct 
stress approximations. However, their stability properties, as a function of both h (the mesh size) 
and p (the polynomial degree of solution approximations) are presently unknown. In this contri- 
bution, we bring to the table, a method of the novel discontinuous Petrov-Galerkin (DPG) type, 
which exhibits stability independent of h and p. The new method is the first /ip-optimal method 
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for linear elasticity that can simultaneously approximate the stress and the displacement. We are 
also able to show, theoretically and practically, that the convergence of the discrete solution does 
not deteriorate as the Poisson ratio approaches 0.5, i.e., the method does not lock. 

To understand how the new DPG method sidesteps the traditional difficulties, let us review the 
usual difficulties in designing schemes that give direct stress approximations. Mixed methods 
[5] based on the Hellinger-Reissner variational principle face the problem of designing an approxi- 
mation space for stresses consisting of matrix functions that are pointwise symmetric. While this 
is not difficult in itself, when combined with two other additional requirements, difficulties arise. 
The first requirement is that forces on a mesh face shared by two mesh elements must be in equi- 
librium. The second requirement is that the conforming stress space, together with a space for 
displacement approximations, form a stable pair. To put this in mathematical terms, let M denote 
the space of real N x N matrices and let § denote its subspace of symmetric matrices. The above 
mentioned stress properties imply that the exact stress a on an elastic body occupying fl C M. N 
lies in the space 

H(div,n-,S) = {ae L 2 (S1;§) : diva e L 2 (f2,R N )}. (f) 

(Here, the set of functions from Q into X whose components are square integrable on Q is denoted 
by L 2 (f2, X), for X = S, M, R N etc.) Mixed methods must use conforming finite element subspaces 
of H(div, i7; S). Although such spaces are known [2], they have too many unknowns (e.g., their 
lowest order space has 162 degrees of freedom on a single tetrahedral element). Such rich spaces 
seem to be necessary to satisfy both the first (conformity) and the second (discrete stability) 
requirement. In contrast, our new DPG methods change the game by separating the approximation 
and the stability issues. 

The DPG method uses a weaker variational formulation for the same problem. In this for- 
mulation, a is sought in the space L 2 ({2,M), in contrast to the space iJ(div, /2;§) above. Since 
L 2 (J?,M), has no interelement continuity constraints, we are able to design an approximating 
finite clement subspace trivially. Furthermore, due to the nonstandard stabilization mechanism of 
the DPG scheme, the discrete stress space can be chosen to be a subspace of L 2 ({2,§), i.e., the 
method gives stresses that are exactly (point-wise) symmetric. One can equally well choose stress 
approximations in a subspace of if(div, fl\ S), disregarding discrete stability considerations. The 
stability of the DPG method is inherited from the well-posedness of the new ultraweak formula- 
tion. Of course, it is by no means trivial to prove this well-posedness (and most of the analysis in 
this paper is devoted to it). It is provable by adopting a Petrov-Galerkin framework where trial 
and test spaces are different. For any given trial space, we can locally obtain a test space that is 
guaranteed to yield stability. 

Test spaces that guarantee stability can be obtained by following the DPG methodology intro- 
duced in [T4"lll5| . Our initial idea was as follows: If one uses DG spaces, then given any test space 
norm, one can locally construct test spaces that yield solutions that are the best approximations 
in a dual norm on the trial space. This dual norm was called the "energy norm" . However, we 
realized |I6j that these energy norms are often complicated to work with once we move beyond 
one-dimensional problems. But we turned the tables in [131134], by showing that given a desirable 
norm in which one wants the DPG solutions to converge, there is a way to calculate the matching 
test space norm. We refer to this norm as the "optimal norm" on the test space (see § 15. 2p . The 
catch is that the optimal norm is nonlocal. Its use would make the computation of a basis for 
the test space too expensive. Hence, we have been in pursuit of norm equivalences. If one uses, 
in place of the optimal norm, an equivalent, but localizable test norm, then the DPG method, 
instead of delivering the very best approximation, delivers a quasioptimal approximation, i.e., the 
discretization error is bounded by a scalar multiple of the best approximation error. This approach 
was applied to a one-dimensional wave propagation problem in |34) . A number of further theoret- 
ical tools were needed to develop an error analysis for multidimensional problems. These, in the 
context of the simple Poisson equation, appear in [13] . In this paper, we further generalize these 
tools to the case of the elasticity problem and introduce new tools to prove locking-free estimates. 

Before we proceed to the details of this DPG method, let us mention several alternative so- 
lutions to handle the difficulty of constructing approximating subspaces of ([lj. One approach is 
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to relax the symmetry constraint of stresses by using a Lagrange multiplier. This means that a is 
sought in 

if(div, J2; M) = {a e L 2 (f2;M) : diva e L 2 (Q,R N )} 

(cf. (HJ), a space for which finite elements are easier to design. This avenue gave rise to mixed 
methods with weakly imposed symmetry [4,9,22,29,30,52]. Yet another avenue is to keep the 
stress symmetry, but relax the i?(div)-conformity. This yielded non-conforming methods, e.g., [3 
23,27!. However, none of these methods have been proved to be hp-optimal. The closest attempt 
to an hp-method is (30) which studies a variable degree mixed method, but does not show how the 
error estimates depend on p. In contrast, we will prove that the (two) DPG methods we present 
in this paper are hp-optimal. Furthermore, since the DPG method can be reinterpreted as a least 
squares method in a nonstandard inner product, it yields matrix systems that are symmetric and 
positive definite (despite having both stress and displacement as unknowns). 

As mentioned above, there are two new DPG formulations in this paper. The first is easy to 
derive and is a natural extension of our work on the Poisson equation in [13] . The second differs 
from the first due to the presence of a scaled Lagrange multiplier. The multiplier serves to obtain 
the extra stability required to prove the locking-free convergence estimates. 

In the next section, we introduce the first DPG method. We also state the main convergence 
result for the first method. The second method and its convergence result is presented in Section|3l 
Then, in Section [4l we study the relationship between these two methods, discovering when they 
are equivalent. The proofs of the above mentioned two convergence theorems appear in Section [5j 
As corollaries to these convergence theorems, we obtain h and p convergence rates in Section [5J 
In Appendix |21 we present a result on a mixed method that we crucially use in our proofs. 



2 The first DPG method 

In this section, we present the derivation of the first of our two DPG methods for linear elasticity. 
We also state a convergence theorem, which will be proved in a later section. 

Linear elasticity is described by two equations. The first is the constitutive equation 

Ao = e(u) (2a) 

and the second is the equilibrium equation 

divcr = /. (2b) 

These equations are imposed on a domain Q C M. N and the space dimension N equals 2 or 3. 
We assume that Q is a bounded open subset of with connected Lipschitz boundary. The 
stress a(x) is a function taking values in S and its divergence (diver) is taken row-wise. The strain 
tensor is denoted by e(u) — (gradu 4- (gradu)')/2 = symgradu where the prime (') denotes 
matrix transpose, and symM = (M + M')/2. The material properties are incorporated through 
the compliance tensor A(x) in (|2a|) which at each x 6 f2, is a fourth order tensor mapping § into 
S. The vector function u : Q h-> M. n denotes the displacement field engendered by the body force 
/ : Q i— > M. N . We consider the simple boundary condition 

u = on df2 (2c) 

which signifies that the elastic body is clamped on the boundary dfl. (We will remark on extending 
the method to other boundary condition later - see Remark [5]) 

To motivate the derivation of the first scheme, we multiply the equations of (2J by test functions 
t : Q \— > S and v : fl ^- > M , supported on a domain K. We temporarily assume r and v to be 
smooth so we can integrate by parts to get 



(A<7,t) k + (u,divr) K - (u,Tn) 1/2jdK = 0, 

{a,Vv) K - (v,crn} 1 /2.dK = {f,v) K . 



(3a) 
(3b) 
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Here n denotes the outward unit normal on the boundary of the domain under consideration, (•, - )k 
denotes the integral over K of an appropriate inner product (Frobenius, dot product, or scalar 
multiplication) of its arguments, and (•, 1)i/2,ok denotes the action of a functional I € H^ 1 ^ 2 {dK). 
Here and throughout, we use standard Sobolev spaces without explanation, e.g., H 1 (f2,M. N ) = 
{v £ L 2 ((2,R N ) : gradv £ L 2 (f2,M)}. 

Now, assume that the domain Q where the boundary value problem ([2]) is posed, admits a 
disjoint partitioning into open "elements" K, i.e., U{K : K £ Qh) — Q. We need not assume that 
elements are of any particular shape, only that the mesh elements K £ flh have Lipschitz and 
piecewise planar boundaries, i.e., in two space dimensions, K is a Lipschitz polygon, and in three 
space dimensions, if is a Lipschitz polyhedron. We now sum up the equations of element by 
element to obtain 

(Aa, r)a h + (u, div r)n h - (it, t ri)dn h = (4a) 
O, Vv)n h ~ (v, a n) dnh = (/, v) f2h (4b) 
(a,q)a h =0. (4c) 

Here, we have additionally imposed the symmetry of the stress tensor by the last equation (where 
q is a skew-symmetric matrix valued test function on ft) and used the following notations: 

(r,s)n h = ^ (r,s) K , (w,l) 9 n h = ^ {w,l) 1/2 ,dK- (5) 
Ken h Kef2 h 

The notation dflh is used for the collection {dK : K £ Qh}- Note that it will be clear from the 
context if differential operators are calculated element by element or globally, e.g., div in (j4|) is 
calculated piecewise, while in (JlJ it is the global. 

The equations of (U) motivate the following rigorous functional framework for an ultraweak 
variational formulation. We let the traces of u and a in the terms (u,t n)QQ h and (v, a tl)qq k to 
be new unknowns, which we call the numerical trace and the numerical flux, resp. The ultraweak 
variational formulation seeks (er, it, u, a n ) in the trial space 

U = L 2 {Q:M) x L 2 {f2; V) x Hg /2 {df2 h ; V) x H-^ 2 {dn h ;Y), (6a) 

satisfying 

b((a,u,u,a n ), (r, v,q)) = 1(t,v, q) V(r,u,g)eF ! (6b) 

where the test space V is defined by 

V = H(div,n h ;§) x H\f2 h ;Y) x L 2 (S7 h ;K), (6c) 

where IK denotes the subspace of M consisting of all skew-symmetric matrices, and the bilinear 
form &(•, •) and the linear form l(-) are motivated by Q. Namely, 

b((a,u,u,a n ), (r,v,q)) = (Aa,r)n h + (u,divr) nh - {u,T<n)dn h (6d) 

+ {(T, ^v)n h - (v, a n ) d Q h + (a, q)n h , 
l(r,v,q) = (f,v). 

In the notations of (|6"a|) and (|6T|) . V denotes the vector space and 

Hl /2 (dU h ;Y) = { V :3w£ H^f2;Y) with r,\ dK = to| a jf,Vif £ O h }, 

H- 1/2 (dn h ;Y)={r)e II H- 1/2 {dK\ V) : 3q £ H (div, f2;M) with 
Ken h 

v\dK = qn\ dK ,\/K £ £2 h }, 
H(div,a h ;S) = {t : t\ k £ H(div,K;§),VK £ Q h }, 
H\f2 h ;W) = {v : v\ K £ H^K-^VK £ f2 h }. 
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The norms on Hl /2 {dQ h ; V) and H- 1 / 2 {dQ h ;N) are denned by 

IHIff 1/2 (d« h ) = ^{INUm") : Vw e #o(^; v ) with « - w \dK = o}, (7) 

\\vn\\ H -m{dn h ) = inf{lkllff(div/2) : Vg G if (div, /2; M) with ay, - gn| 9A ' = 0}. (8) 
The trial and test norms are defined by 

\\(a,u,u,a n )\\^ = Ho + Nlo + Nl^/'^) + \\^n\\ 2 H -i/2 {dnh) , (9a) 

\\(r,v,q)\\ 2 v = ||r[& (diVilJh ) + \\v\\ 2 HHnh) + \\q\\ 2 n . (9b) 

Here the norms on if 1 (J7/ t ; V) and H(div, S?h,;S) are defined by 

\\ v \\m(n h ) = i v , v )n h + (gradv, gradv)o h , (10) 

IMIff(drv,f2 h ) = {T,T)n h + (divr,divr) r2fe . (11) 



This completes the description of our new infinite dimensional DPG variational formulation. 

Remark 1 Like in many other numerical formulations, in the DPG formulation ([5]), we need to the 
extend the domain of the compliance tensor A(x) from S to M. There are many ways to perform 
this extension. To choose one, decompose M orthogonally (in the Frobenius inner product) into K 
and S. A standard way to extend A(x) from § to M is to define A(x)n — k for all k in K. Then, 
whenever the original A(x) is self-adjoint and positive definite on §, the extended A(x) is also 
self-adjoint and positive definite on M. 

We assume throughout that A(x) (i.e., its above mentioned extension) is self-adjoint and pos- 
itive definite uniformly on Q. We also assume that the components of (the extended) A are in 
L°°(f2). 

Next, we describe the discrete DPG scheme. This is done following verbatim the abstract setup 
in [T31 Section 2] (see also [H]). Accordingly, we define the trial-to-test operator T : U n> V by 

{TU,V) V =b(U,V), W G V and \JU G U. (12) 

We select any finite dimensional subspace Uh C U and set the corresponding finite dimensional 
test space by 

V h = T{U h ). 

Then the DPG approximation (crf l ,U} l ,Uf l ,d'n,h) G Uh satisfies 

b((ah,u h ,Uh,d- n!h ), (r,v,q)) = l(r,v,q) V (r, v, q) £ V h . (13) 

The distance between this approximation and the exact solution can be bounded as stated in the 
next theorem. 

Theorem 1 (Quasioptimality) Let Uh C U. Then, (I6bp has a unique solution (a,u,u,a n ) G U 
and (|13p has a unique solution {o~h,Uh,Uh,a nt h) & Uh- Moreover, there is a > independent 
of the subspace Uh and the partition flh such that 

V < C (1) A, 

where T> is the discretization error and A is the error in best approximation by Uh, defined by 

T>=\\a- o- h \\ L 2 {n) + \\u - u h \\ L 2 {n) + \\u - u h \\ H i/2 {dnh) + ||<7„ - CTn,h||H-i/a(an h )- 
A = inf 

Mk - Ph\\mn) + \\u - w h \\ L 2 {n) + \\u - Zh\\ H i/2 {dnh) + \\a n - Vn,h\\H-y^(dQ h ) 
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This result is comparable to Cea's lemma in traditional finite element theory. Of importance 
is the independence of with respect to Uh- Specifically, we are interested in setting Uh to 
/ip-finite element subspaces with extreme variations in h and p to capture singularities or thin 
layers in solutions. In this case, the constant being independent of Uh, is independent of 

both the mesh size h and the polynomial degree p. As such, this forms the first method for linear 
elasticity with provably hp-optimal convergence rates of the same order for a and u. Although 
several mixed methods yielding good approximations to a are known, proving their /ip-optimality 
requires proving an inf-sup condition carefully tracking the dependence of constants on p, a feat 
yet to be achieved. For a proof of Theorem [TJ see Section [5] 

3 The second DPG method 

The robustness of numerical methods with respect to the Poisson ratio is an important consider- 
ation in computational mechanics. Methods that are not robust exhibit locking. Note that we did 
not assert in Theorem Q] that the constant C^ 1 ' is independent of the Poisson ratio. However, in 
all our numerical experiments (see Section [7} , the method showed locking- free convergence. This 
section serves as a first step towards explaining this locking-free behavior theoretically. 

The second DPG method given below is designed so that we can establish locking-free con- 
vergence with respect to the Poisson ratio. It has one more trial variable. In this section, we will 
assert its locking-free convergence properties, restricting ourselves to isotropic materials. In the 
next section, we will provide a sufficient condition under which the first and the second DPG 
methods are equivalent. This gives theoretical insight into the locking-free behavior of both the 
first and the second methods. 

Let us begin by defining the essential infimum 

Qo = essinf (tr(A(x)I)). (14) 

Obviously, Qq > for an isotropic material with Poisson ratio v < 0.5. The second method is 
motivated by the same integration by parts as in Q , but with the following additional observation 
in mind: If we set r — I in (U) and recall that u = Q\on, then we have that 

/ tr/ler = 0. (15) 
J n 

Imposing this condition via a Lagrange multiplier, we obtain another ultraweak formulation with 
the following bilinear and linear forms: 

b {2 \{cr,u,u,& n ,a), (r,v,q,f3)) = (Ao-,r)n h + (u, divr)f 2(i 

+ (a, Vv)n h + (a,q)n h 
+ Qa 1 {Aa,pI)a h 
l^(T,v, q ,f3) = (f,v). 

Here / is the identity matrix. The trial space is now set to 

U (2) = L 2 (tt;M) x L 2 (J?;V) x Hg /2 (df2 h ;W) x H- 1/2 {dQ h ;Y) x R (16b) 

and the test space is set to 

y( 2 ) = ff(div,/4;S) x H x {fi h \Y) x L 2 (f2 h ;K) x K. (16c) 

By (fT5"l) , the solution (a, u) of will be the solution of (fT5|) with 

u = u\ d n hl v n = o-n\dQ h , and a = 0. 

However, more work is needed to conclude similar statements at the discrete level (see the next 
section). 



{u,Tn) anh + Q 1 (aI,AT) f2h 
- (v, o- n )dn h 



(16a) 
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Remark 2 In the case of mixed boundary conditions we must add the term 

-Qo 1 (fi,C8i)n)sn h (17) 



to the the expression (|16a[) . Note that this term vanishes in the case of kinematic boundary- 
conditions analyzed in this paper. However, for more general boundary conditions, u can be be 
nonzero on the parts of the boundary where traction conditions are imposed. Hence (|17|) simplifies 
to a boundary integral that is nonzero in general. 

The second DPG method is obtained by constructing a discrete scheme as before from the 
ultraweak formulation (following [13l Section 2]). The trial-to-test operator in this case (cf. (flBl ) 
is T ( 2 ) : i — y VW by 

(T^U,V) V (2) = b (2 \u,V), VVey (2) . (18) 
Let Uj® C U {2) be any finite dimensional subspace. We set = T^iU^). The second DPG 

(2) 

approximation [ah,Uh,Uh,o'n,h> a h) S U h satisfies 

bW((a h ,u h ,u h ,a nih ,a h ), (r,v,q,P)) = l((r,v,q,P)) V (r, v, q, (3) G vf } . (19) 

As in the case of the first DPG method, we are able to prove a quasioptimality result (see the 
next theorem) bounding the discretization error 

£> (2) = ||<7 - (7 h \\ L 2 {n ) + \\u - u h \\ L 2 {n) + \a- a h \ 

+ l|w - Uh\\ H V2 {dnh) + \\a n - a nth ||fl--i/2(ao h )- 

by the error in best approximation 

A {2) = inf I ||cr - phWwn) + \\ u - w h \\ L 2 {n) + \a - j h \ 

(Ph,Wh,Zh,fln,h,lh)£U h \ 

+ II" - *h\\ H i/2 {gs - 2h) + \\a n - f}n t h\\ H -y»(dn h ) 

However, we are also able to prove a stronger result under the following assumption. 
Assumption 31 (Isotropic material) We assume that 

At = Pt d + Q^I (20) 



where 



tr(r) 



for any r in M, and P and Q are positive scalar functions on fi. (Then, obviously Q > Qq for 
the Qo defined in (|14[) . ) When (|20[) holds, we also define 

B = Q^\\Q\\ Lx(n) , (21) 
P = ess inf P(x). (22) 

(Note that (|20[) is assumed to hold for all r 6 M..) 

Remark 3 When we consider isotropic materials, we do not extend A from § to M in the way 
suggested in Remark [1] Instead, we assume that it is extended from S to M by An = Pk for all 
K6l. This ensures that {2D) holds for all r e M, not just for all r e §. 
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Theorem 2 (Quasioptimality of the second DPG method) Let U ( h 2) C . Then, (H]) has 
a unique solution (a,u,u,a n ,a) £ £/( 2 ) and (fTO]) has a unique solution (ah, Uh,Uh,a n ,h,cth) <E . 
Moreover, there is a C*- 2 - 1 > independent of the subspace U? and the partition flh such that 

£>(2) < (7(2) ^(2) _ 

If in addition, Assumvtion \31\ holds, then can be chosen to be 

= cP \\\A\\ + BfB\\\A\\ + P + l) 2 , (23) 

a constant independent of Qo, and consequently the method does not lock. (Here, the positive 
constant c is independent of A.) 

The proof of this theorem appears in Section [5l 
4 The relationship between the two DPG methods 

In this section we will establish that for homogeneous isotropic materials the two DPG methods are 
equivalent. We will also show that despite the additional trial variable, the second DPG method 
can be solved in essentially the same cost as the first. 

4.1 The equivalence 

Recall that the first DPG variational formulation uses the trial space 

U = L 2 (Q;M) xL 2 (f2;Y) x H X/2 (dn h ;Y) X H- l / 2 (d{2 h ; V), 
while the second uses U x K. We begin with the following simple lemma. 
Lemma 1 For any U = (a,u,u,a n ) in U, let 

(r,v,q) = TU. 

Then, with = (a,u,u,a n ,0), 

T^U^ = (r,v,q,P), 

with $ = Qo X {Aa, I) Q . 

Proof By the definition of T, we have, for any (S T ,S v ,5 q ) E V = H(div, fl h ; §) x H 1 (fi h ;N) x 
L 2 (f2 h ;K), 

( (r,v,q), (S r ,S v ,Sg)) v = (t,5 t )o + (div r, 5 T )n h + (v,6 v )q + (Vv,V5 v )a h + (q,5 g )n 

= (Aa, 5 T )n + (u,divS T )n h - (u,S T n) d n h + (o-,V5 v )a h - (S v ,a n ) d n h + (a,8 q ) n . 

Therefore choosing 

(i^Q^(Aa,I) n , 

we obviously obtain 

( (r,v,q,0), (6 T ,5 v ,6 q ,6p)) v(2) 

= {t,S t )(2 + (divT,8 T )n h + (v,S v )n + (Vv,W6 v )a h + (q,S g )n + 

= (Aa, 5 T )n + (u, div S T )n h - (u,S T n) d a h + (a,X75 v )n h - (S v ,a- n }dQ h 

+ (tr, S g )n + Qa 1 (Aa, 8 t3 I)n 

= b [2) ((a,u,u,a n ,Q), (r,v,q,j3)) 

for any (S T ,S V , 8 q , Sp) G V x K = V^ 2 \ This finishes the proof. 
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We use the above result together with the following assumption to prove the equivalence. The 
assumption essentially states that the material is homogeneous and isotropic and that the discrete 
trial space at least contains two specific functions related to the identity matrix. 

Assumption 41 Suppose Assumvtion \31\ (on isotropy) holds. Let the discrete trial subspace Uh C 
U of the first DPG method be 

U h = E h xW h x M h x F h , (24) 
and let the second DPG method use the trial space Uh X K. We assume that 

1. Q(x) — Qq for all x G CI, 

2. I G S h , 

3- In\ 9 Q h G F h . 

Lemma 2 If Assumption\4l\ holds, then (7,0,0) G Vh = T(Uh)- 

Proof By virtue of the assumption, the trial function U — (o~,u,u,a n ) with a = I, u = 0, a n = 
In\ga h , and u = 0, is in Uh- Hence, (r,v,q) = TIA is in Vh and satisfies 

[t,8 t )q + (divT,div(5 T )f 2ft + (v,5 v )q + (Vv,yS v ) i2h + (q,5 q )n . . 

= (AI, 5 T )n + {I, VS v )n h - (S v ,In) 9 n h + (I, S q ) n 

for all (5 T ,5 V) 8q) eV = H(div,n h ;S) x ff 1 ^; V) x L 2 (S7 h ;K). The last term in ^ vanishes 
due to the skew symmetry of S q . Integration by parts shows that (I,S75 v )s-2 h — (5 v ,In)QQ h = 
also. Hence, we conclude that (r, v, q) = [QqI, 0, 0) is the unique solution of (1231) . Since (t, v, q) is 
in Vh, we have proved the lemma. 

Lemma 3 // Assumption \Jl\ holds, then 

(7,0,0,0) =T( 2 '(0, 0,0, 0,1). 

Proof Let (r, v, q, (3) = T< 2 )(0, 0, 0, 0, a). Then by the definition of , we have 

(t, S T ) n + (divr, div(5 r ) f 2 h + (v,6 v )q + (Vv,V5 v )i 2h + (q,5q)n + PSp 
= Qo 1 (aI,AS T ) n 

for all (5 T ,S v ,S q ,Sp) G V^ = H{div, f2 h ;§) x il^fi^V) x L 2 (!7 h ;K) x M. Putting a = 1, and 
using the symmetry of A, the right hand side Q^ 1 {al, A5 t )q = Qq 1 (AI ,8 t )q, which due to the 
assumption on A equals (J, 8 t )q, i.e., 

(t,S t )q + (divT,divS T )n h = {I,8 t )q. 

It is now obvious that (r, v, q, (3) = (I, 0, 0, 0). 

Theorem 3 Suppose Assumption \4-l\ holds. Then (o~h, Uh, Uh, o ni h) solves the first DPG method (113|) 
if and only if (ah, Uh, iih, o n ,h, 0) is the discrete solution of the second DPG method (fTS)) . 

Proof Suppose (o~h,Uh,Uh, <r n ^h) solves the first DPG method ([6]). To show that the function 
(<j/i, Uh, Uh, &n,h, 0) satisfies the second DPG method (|16l) . we have to show that the equations 



(Ao h ,T)n h + (uh,divT)n h - (u,Tn) di2h + Q^ial , Ar)n h = (26a) 

(o- h , Vv)n h - (v, a n ,h)dn h = (/, v) (26b) 

(a h ,q)o h =0 (26c) 

Qo\Aa,f3I) nh =0, (26d) 



hold, with a — 0, for all (r, v, q, j3) G T^(Uh xl). To this end, we proceed in two steps. 
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First, we consider test functions of the type T^ 2 '(Uh X {0}). By virtue of LemmaQ] these test 
functions are in T(Uh) x R. Hence by the fact that (ah, Uft, Uh, &n,h) solves the equation of the first 
DPG method (USD for all (r,v,q) £ T(U h ), we observe that (ggg) , (j26b|) and pjcl) hold. To show 
that (I26dj) also holds, we observe that because of Lemma [3] we may put (r, v, q) = (/31, 0, 0), for 
any /3 S R, in (13]) to get 

(Aa h ,j3I)n h + (u h ,div (3I)n h - (u h ,filn) d n h = 0. 

Multiplying by Qq 1 and simplifying, we obtain (|26d|) for any fi. 

Second, consider test functions of the type T^ 2 ^({(0, 0, 0, 0)} x R). But by Lemma [3] we know 
that T( 2 )(0, 0,0,0, 1) = (7,0,0,0). So to show that (gBJ) holds for this test function, it suffices to 
prove that (|26a|) holds with t = I. But this follows from LemmadJ which shows that (I, 0, 0) is in 
V h . 

Conversely, suppose (<Jh,Uh, Uh, an,ht 0) satisfies the second DPG method (|19p . Then, (ah, Uh, Uh,o n ,h) 
must satisfy the first DPG method, for if not, there must be another function (a' h , u' h , u' h , a' n h ) solv- 
ing the first DPG method. But then, the already proved implication shows that (a' h , u' h , u' h , a' n h , 0) 
must solve the second DPG method. This contradicts the unique solvability of the second DPG 
method asserted in Theorem [2] 



4.2 Solving the second DPG system 

The relationship between the two methods revealed above, can be utilized to obtain a useful 
strategy for solving the second DPG system. 

To understand the linear systems that result from both methods, we assume that a basis for 
Uh is made of local (standard finite element) functions 

e-i = (fi, Ui, iii, a n ,i) 

for i — 1,2, ... ,m. Then the corresponding basis for Vh is furnished by tj = Te{. Hence the m x m 
stiffness matrix of the first DPG method E has entries 

Eij = b(ej,ti). 

It is easily seen that this matrix is symmetric and positive definite. (This is a general property of 
DPG stiffness matrices - see [13] or [15].) In addition, E is sparse due to the locality of the basis 
functions. 

Now, consider the stiffness matrix of the second DPG method. Here, we need a basis for Uh X R. 
It is natural to take as basis for Uh x R, the m + 1 functions defined by 




(ei, 0), for all i = 1,2, ... , m, 

(0,0,0,0,1), fori = m + l. 



But we should then note that the stiffness matrix 

E^=b^(ef,T^ef ) ) 

is no longer sparse. This is because T^ef^ is not locally supported. Indeed, if we write 

T (2) (2) 

as (Ti,Vi,qi, Pi), then can be globally supported even if ef^ is local. Thus, although E^ is 
symmetric and positive definite, one may be led into concluding that the second DPG method is 
too expensive due to the non-sparsity. 

However, this is not the case. Below, we will show how to solve a system E^x^ = j/ 2 ' by 
solving a system Ex = y and performing a few additional inexpensive steps. A key observation is 
that E is a rank-one perturbation of E. 
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Proposition 1 Decompose the matrix E^ into 

= 

where c G W n . Then 



E c 
d d 



E = E + M 
where £ G R m is defined by tj = Qg 1 (Aa jl I) n . 

Proof Let (ri,Vi,qi) = Te^. Then, by Lemma[TJ T^ef^ = (Ti,Vi,qi,/3i) with 

ft =Qo 1 (Aa i ,I)o- 

Hence, 



E ij = b^{ef\T^ef ) ) 

= 6« (ef , (n, v i ,q i ,0))+ (ef\ (0, 0, 0, ft) ) 

= b(e j ,Te i ) + Qo 1 (P i I,Ao- j ) n 
= E ij +p i p j . 

The result follows because ft = li. 

A consequence of this proposition is that we can invert E using the Sherman-Morrison for- 
mula [55], namely 



(E + ££')- = E- L - a(E- i e)(E- L £y, with a 



1 



1 + i'E-H' 



(27) 



Therefore, to conclude this discussion, consider the matrix system arising from the second DPG 
method (fHfl): 



~E c 




X 




fJ 


d d 




JJ_ 








Since E is nonsingular, the solution is given by 

y = -{d - dE- 1 c)- 1 dE- 1 g, x = E~ l g - E^cy. 

Hence the solution of the second DPG method can be obtained by solving the two linear systems 
Ex c = c and Ex g = g. Each of these systems can be solved using formula (|27p at essentially the 
same cost as solving a system involving E. 



5 Error analysis 

In this section, we prove Theorems [T] and [2J We will give the proof in full detail for the more 
difficult case of Theorem [5] first. Afterward, we will indicate the minor modification required to 
prove Theorem Q] in a similar fashion. The plan is to use the abstract DPG framework developed 
in [T3lH4l[T5lfT6 , 34 . summarized next. 
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5.1 Abstract quasioptimality 

Let X with norm || • \\x be a reflexive Banach space over R, Y with norm || • ||y be a Hilbert 
space over R with inner product {-,-)y, and &(•, •) : I x y h> M be a bilinear form. The abstract 
trial-to-test operator T : X M> Y" is defined - as before - by (Tu,v)y — b(u,v) for all u G X and 
y € Y. We can write the DPG method using an arbitrary subspacc Xh C X even in this generality. 
Let y h = T(X h ). 

Theorem 4 (see [13, Theorem 2.1]) Suppose x G X and Xh G Xh satisfy 

b(x,y)=l(y) VveY, 
b{x h ,y) = l(y) VyeYfc. 

.Assume </iai 

{w G X : 60, u) = 0, Vw G F} = {0} (28) 
and i/iat there are positive constants C\ , C2 such that 

Ci\\y\\ Y <\\y\\ optx <C 2 \\y\\ Y , VyeY (29) 

where the so-called "optimal norm" is defined by 

11 11 b(x,y) 

IMI op t,Y = SU P. \(ZnZ ■ ( 30 ) 

Then 



opt ' Y xwex \\x\\x 



C 2 

Xh\\x < jr inf ||a;-Zfc||x. 



5.2 The optimal test norms of both methods 

We now see what is the norm defined by (|3U|) in the context of the first and second DPG methods. 
From the structure of the bilinear form it is easy to see that, for the first DPG method, the optimal 
test norm is 

b((o-,u,u,a n ), (T,v,q)) 2 

0yi(a,u,u,a n )£U || (tT, U, M, <T„) ||^ 

= ||Ar + Vi; + a|| 2 r2h + ||divr||^ + ||[rn] 
where the "jump" terms are defined by 



( T ,v,q)\\ 2 opty = sup 



(31) 



\d(2 h 



\[vn]\\aa h ■= 



(u,Tn)on h 

su p in 



sup 



(v, a n) 



an h 



(32a) 
(32b) 



0#erGff(div,fi;M) IMI.ff(div,fi) 

Similarly, for the second DPG method, taking the supremum over its trial space, we have 

\\(T,v,q,P)\\l P t,vw =\\AT + Vv + q+Qv l pAI\\l h + \\&VT\\l h (33) 



+ rn 



\ao h 



\\[vn] 



\ao h 



[ tr(Ar) 
Jn 



In either case, the optimal norms are inconvenient for practical computations, due to the last two 
jump terms. These terms would make the trial-to-test computation in (jl8p non-local. Therefore, 
a fundamental ingredient in our ensuing analysis is the proof of equivalence of the optimal norm 
with the simpler "broken" norm (for the first DPG method) 



II (r,v,q) f v = ||r|& (dlViflfc) + \\vf m{ a h) + \\qf„ 



(34) 



A DPG method for elasticity 



13 



which does contain the jump terms. For the second DPG method, we will similarly prove that (|33p 
is equivalent to following norm 

||(r,«,g^)||^ w = \\T\\% {d ^ nh) + \\v\\ 2 HHnh) + \\q\\ 2 n + \p\ 2 . (35) 

These equivalences would verify condition (|29p . so we would be in a position to apply Theorem 2J 
Before we proceed to prove these, let us verify the other condition ([25)) . We begin with a necessary 
preliminary. 



5.3 Korn inequalities 

We will need two well-known inequalities due to Korn. The first Korn inequality asserts the 
existence of a constant C > such that 

IM|j*i(fl)<C||e(i;)||n to effort), (36a) 
while the second Korn inequality gives a constant C > such that 

IMUi(«) <C{\\v\\ n + \\e{v)\\ Q ) VveHW. (36b) 

Above and in the remainder, we use C to denote a generic constant whose value, although possibly 
different at different occurrences, will remain independent of the discrete approximation spaces. 
These inequalities can be found in many references. E.g., for (|36a[) . see [26l eq. (2.7)], and for (|36b|) . 
see [55J Section 1.12 in Chapter 6]. 



5.4 Uniqueness for the second DPG method 

In this subsection, we verify condition (|28p for the second DPG method. 

Lemma 4 With and as set in (jTSJ) suppose (a,u,u,a n ,a) G satisfies 

bW((a,u,u,a n ,a),(T,v,q,F3))=0, V(r, v, q, 0) G V< 2 >. (37) 

Then {a, u, u, a n ,a) = 0. 

Proof Equation (|37p is the same as 

(Aa,T) K + {u,divT) K -{u,Tn) aK + Qo 1 (aAI,T) K =0 Vt g H (div, K ; S) (38a) 

{a,gradv) K -(a n ,v} dK = Vu G H x (K; V) (38b) 

tr(Ar) = 0. (38c) 



Here, cr G L 2 (^;§) because 

(cr, g ) r2 ^0 VgeL 2 (^;K). 

We take r G T>(K;§) and u G T>(K;V) arbitrarily. Here, as usual, T>(D, X) denotes the space of 
infinitely smooth functions from D into X that are compactly supported in D. Then we have 

Aa - e(u) + Q 1 aAI = in K (39a) 
div a = in if (39b) 

in the sense of distributions, for every K G fit- In particular, this implies that u G 7? 1 (_ftT;V) by 
f (|39ap and the second Korn inequality (I36bp 'l and a G H (div, if; S) by (I39b|) . 
Now we claim that we also have 

a n \dK = o-n\ dK and u\ dK = u\dK- (40) 
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The first identity in (j4"0)l is obtained by integrating (|38b[) by parts and using (|39bp to find 
that (v,a n — o"n) 1 / 2: QK = for every v 6 i? 1 (J7/ l ;V). Note that this identity implies that 
a £ H(div,f2;S). 

To prove the second identity of (f4"0)) . we need to proceed a bit differently. Let F be an arbitrary 
face of K, say in three dimensions (the two dimensional case is similar and simpler). We will now 
show that given any rj = £ V(F;Y), there is a r 6 H(div, K;E>) such that rn £ L 2 (K;V) is 
supported only onFCK such that tti\f — r\. To this end, we may, without loss of generality, 
assume that F is contained in the xy-plane, so that n = (0, 0, 1)'. Then, set 



r]i{x,y) 

o o m(x,y) 
vi(x,y) m{x,y) m{x,y) 



(41) 



where \ is an infinitely smooth cut-off function such that (i) the support of x(ie, y, 0) is compactly 
contained in F and contains the support of 77, and (ii) \ = 1 on the support of 77. Clearly we can 
find such a cut-off function, and furthermore, construct it so that rn vanishes on all other faces 
of K. 

We use such r to prove the second one in (|40|) . First observe that from (|57|) . we have (Aa, t)k + 
(w, div t)k — (u : t u)qk + Qq 1 (cxAI, t)k — 0, for all r e iJ(div, K; S). Integrating by parts (which 
is permissible since by (1551) . u e H l (K;Y) and a 6 iJ(div, K; §)), and using (|38ap . 



(tt-u,rn) K = 0, Vr e i7 (div, if ; §). (42) 
Choosing r as in (|4"Tj) . this implies 

(u - u) ry ds = 0, Vr? G V(F,Y). 

F 

Hence = in L 2 {F) and this holds for all faces of K. This proves the second identity of (|4U)l . 
which implies that u £ Hq(S1;V) (after also noting that w|a,Q = 0). 

Next, we choose r = 7 on Q and sum the terms in (|38a[) over all if £ .!?/,,. Using the fact that 
u\dK = u\qk, we have that 



/ tr(Ar) + Qq l a [ tr(AT) = 0. 



The first term vanishes due to p8cp . Hence we have shown that a = 0. 

Now, choose r = a and u = u in (|38a[) - (l38bp . Summing up these equations and canceling terms 
after integrating by parts, we find that 

(Aa, <j)n h + (u, a n) dnh - (it, a n) d n h - (u, v n )dn h = 0. (43) 

Since u £ iig(!7;V) the last term on the left had side vanishes. Furthermore, since we already 
showed that the interelement jumps of an are zero, the penultimate term (u, a n)QQ h also vanishes. 
Note finally that (u,an)dn h — (u,an)Qf2 — as u £ H^(Q). Thus, (14U1) implies that a = 0. 

Since a vanishes, by (|39ap . we have that e(u) — 0. Since u £ Hq(Q;Y), by the first Korn 
inequality (|36aD we find that u = 0. Since both a and u vanish, by (f4"0|) . u and ct„ also vanish. 



5.5 An inf-sup condition 

In this subsection, we verify that the lower bound in the condition (f2"9"]l holds for the second 
DPG method. The lower bound is the same as an inf-sup condition due to the definition of the 
optimal norm. To prove this inf-sup condition, we use a modification of the mixed method for 
linear elasticity with weakly imposed symmetry, given in Appendix [21 
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Lemma 5 There is a positive constant C\ such that for any (r,v,q,/3) £ V^ 2 \ 

Ci\\(T,v,q,(3)\\ V (2) < \\{T,v,q,f3)\\ opt v{2) . 
If, in addition, Assumvtion \Sl\ holds, then the constant C\ can be chosen to be 

Cf 1 = c 2 P r 7 1 (||A|| + BfB\\\A\\ +P + I) 2 , (44) 
where c-2 is a positive constant independent of A. 

Proof Given (r, v, q, (3) G V^ 2 \ we solve the mixed method ([70)1 with data F\ = r, F2 = —V, F3 = q 
and F 4 = N-^I/l^l, to get (a,u,p,a) G iJ(div, Q\ M) x L 2 (J?;V) x L 2 (f2;K) x E. According to 
Theorem [51 the component u is in Hq(Q;Y) and 

IMIff(div,n) + IIwIIh 1 ^) + + M < c ^\\ T \\n + + \\q\\a + |/3|). (45) 



Within an element K, we can find the strong form of the equations in (1701) by choosing infinitely 
smooth test functions that are compactly supported on K. We obtain 

Acr - Vit + p + aQo 1 AI = r, (46a) 

div a — —v, (46b) 

skw a = q, (46c) 

together with the last equation (I70d[) which can be restated simply as 

Qo 1 / tr(Ao-) = /3. (46d) 



These, together with integration by parts, imply that 

\\r\\ 2 n + \\v\\ 2 n + \\q\\l + \p\ 2 

= (Acr - V« + p + aQ^AI, t)q - (div a, v) n + (skw a, q) n + (3Qq 1 / tr(Aa) 

Jn 

= (a, Ar)n h + {u, divr)^ - (u, rn} df2h + {p, r) n + aQ 1 / tr(Ar) 

J n 

+ (a,Vv)a h - (v,an) d Q h + (skwa,q) n + (cr, PQq 1 AI)a h 

= (a,Ar)a h + (u, divr)^ - (u,Tn) d Q h + ciQ ~ 1 / tr(Ar) 

Jn 

+ (cr, Vv) i2h - (v, o-n) 9 n h + (c, q)n + {<r, /3Qg 1 AI)n h 
since r is symmetric. Rearranging and applying Cauchy-Schwarz inequality, 

\\r\\l + \\v\\l + \\q\\l + W\ 2 

= {a,AT + \7v + q + l3Q^AI) nh + (u, div T ) Qh + aQ^ 1 / tr(Ar) 

Jn 

- {u,rn) dnh - {v,an) d a h 

< \\o-\\ n \\AT + Vv + q + pQ^AlWn, + \\u\\ n \\ divr|| fih + \a\ ■ [ tr(Ar)| 

Jn 

+ \\[Tn}\\ ds - 2h \\u\\ H i in) + \\[vn]\\ d n h \\a\\ H(diVt0) 

< 2 \\(T,v,q,(3)\\o P t,vW (IMl!r(div,n) + \H\ 2 H1{n) + \\p\\% + \a\ 2 ) 1/2 
By (gS]), we have that 

\\r\\% + \\v\\l + \\q\\ 2 n + \(3\ 2 < 4C 2 \\(r, v,q,(3)f op t,vW ■ (47) 
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Furthermore, since \\At + Vv + q + f3Q 1 AI\\a h < ||(r, v, q, (3) || opt v(2) , by triangle inequality, 

||V«|| flh < ||A||||r[[ n + \\q\\a + IIATHQo 1 ^] + \\(r 7 v 7 q, (3)\\ opt ,v^ ■ 
which implies that 

[[Vu||fl„ < c-t\\(T,v,q,P)\\ opty w. 

for a positive constant c\ . 

If, in addition the material is isotropic in the sense of Assumption [2U then using the notations 
of the assumption, the constant c\ can be chosen to be 

ct = 2(|L4|| +B)C . 

Finally, since || divr||f2 h < ||(t, v, q, ft) || t y<2), we can control all terms in forming the norm 
||divr|| fih , i.e., 

\\{T,v,q,p)\\ v m < c 2 (2C +c 1 )\\(T,v,q,l3)\\ optyi 2) 

with a constant c 2 is independent of A. The lemma follows with Cj" 1 = C2(2Co + c\). In the case 
of isotropic material, observe that 

Cr 1 = c 2 (2C + c x ) = 2c 2 (||A|| + B + 1)C < 4c 2 (p|| + B)C 
< 4c 2 (\\A\\ + B)c 1 P - 1 J B 4 (||^|| + P + 1) 2 (||A|| + B) 
Kc^^dlAH+B^dlAII+Po + l) 2 , 

where c 2 is a positive constant independent of A. 
5.6 Upper bound 

Now we show that the upper inequality of condition (|2"5j) can be verified for the second DPG 
method. 

Lemma 6 There is a positive constant C 2 such that for any (r,v,q,(3) £ V^ 2 \ 

Ufa v > 0) Hopt^W < c 3 (\\A\\ + B)\\(t, v, q, p)\\ vw . 
Here, C3 is a positive constant independent of A. 

Proof Let us first prove an upper bound for the jump terms. Integrating by parts locally and 
applying Cauchy-Schwarz inequality, 

llr ill (w,Tn) dnh (gT&dw,T)n h + (w,divr) fih 

[ rn ] \\dn h = sup — n — n = sup n — n 

" weH^(n-,v) \M\Hi(n) weH^(n-v) \\w\\h^{0) 
< \\ T \\H(di v ,n h )- 
We use a similar argument for the other jump, i.e., 

||[ H || = sup pSHkB, = Bup (gradt,0„ h +KdivO nh 

h ?eff(div,J?;M) IKIIff(div,i?;M) <;eff (div : l?;M) 11? II ff(div,fi;M) 

< \\v\\m(n h )- 
The remainder of the proof is straightforward. 



5.7 Proof of Theorem [T] 



We apply the abstract result of Theorem |U Assumption (12"51) is verified by Lemma [U The lower 
inequality of (|29p is verified by Lemma [5j and the upper inequality is verified by Lemma El Q 
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5.8 Proof of Theorem H 

The analysis of the first DPG method is in many ways simpler than the above detailed analy- 
sis of the second DPG method. Just as in the proof of Theorem [TJ we only need to verify the 
conditions (|28|) and ([29)) of the abstract result. 

The proof of the uniqueness condition (|28p is similar and simpler than the proof of Lemma 21 
so we omit it. 

The proof of the upper inequality in condition (1291) is the same the proof of Lemma [51 
The proof of the lower inequality in condition (|2"9"]l for the first DPG method is analogous and 
simpler than the proof of Lemma [5j To highlight the main difference, instead of considering the 
mixed method in Appendix we now need only use the standard mixed method with weakly 
imposed stress symmetry. In other words, the analogue of (|45[) is now obtained as follows: Given 
(t, v, q) G V, we solve the following variational problem to find (a, u, p) G H (div, Q\ M) xL 2 (i7; V) x 



L 2 (i?;K) satisfying 

(Aa, St) n + (u, div St) n + (p, St) q = (t,St) q VSt G ff(div, !?;M), (48a) 

{diva,Sv)o = -(v,Sv) n VSv E L 2 {Q;V), (48b) 

(<T,Sq)n = (q,Sq) n VSq E L 2 (H;K). (48c) 

Then we use the standard stability estimate [3] for this method to get the analogue of (|4"5")l and 
proceed as in the proof of Lemma [5l □ 



6 Examples of trial spaces and convergence rates 

The trial subspaces of both the first and the second DPG methods (namely Uh and U h ) were 
unspecified in Theorem [T] and theorem [51 This section is devoted to two examples of trial spaces 
and how one can use Theorems [T] and [2] to predict h and p convergence rates for these examples. 
The examples we have in mind are DG spaces built on a tetrahedral mesh and a cubic mesh. We 
only consider the case of the first DPG method (as the same convergence rates can be derived 
analogously for the second DPG method). 

If D is a simplex, let P P {D) denote the set of functions that are restrictions of (multivariate) 
polynomials of degree at most p on a domain D. If I? is cubic, then we write it as a tensor product 
of three intervals D = D x ®D y ®D z and define Q p ' q ' r (D) = P p (D x )®P q (D y )®P r (D z ) which is the 
space of polynomials of degree at most p, g, r with respect to x, y, z, resp. As with Sobolev spaces, 
when these notations may also be augmented with a range vector space, i.e., P p (D;E>) denotes 
the space of symmetric matrix-valued functions whose components are polynomials of degree at 
most p, etc. 

Recall that - sec (124)) to specify Uh, wc must specify its four component spaces. If Qh is a 
tetrahedral mesh, we set 

£h,p = {p:p\K£P P (K;§)}, W h>p = {v : v\ K e P P {K; V)}, (49a) 

while if flh is a cubic mesh, then we set 

Z h , P = {p: P \k eQ p ^ p (K;§)}, W h , p = {v : v\ K £ Q P > P ' P (K;Y)}. (49b) 

Note that we have chosen the subspaces to consist of symmetric matrix polynomials. This is 
clearly allowed since the only requirement for the discrete trial subspace was that Uh C U. (The 
corresponding automatically generated test space ensures stability of the resulting method. We 
emphasize that this stabilization mechanism is different from mixed methods.) The numerical flux 
space is set as follows: 

Fh, P — {v '■ v\e € P P (E;Y) V mesh faces E} if flh is a tetrahedral mesh, (49c) 

F h , P = {q: v\e G Q p ' p (E; V) V mesh faces E} if Q h is a cubic mesh. (49d) 
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In either case, we define the numerical trace space by 

M h . p+ i = {r]:3w £ W htP+ i n H^(Q;Y) such that tj\ k = w\ K VK E Q h }. (49e) 

Since p > 0, the space Mh, p +\ is non-trivial. 

Let us apply Theorem [2] with these as trial spaces for each solution component. Then, if we 
know how the best approximation error converges in terms of h and p, we can conclude rates of 
convergence. It is well known that for s > 0, 

inf \\u-w h \\ n < Ch s p^ s \u\ H s (n) , (s<p + l). (50) 

Here p 2 = max(p, 2). A similar best approximation estimate obviously holds for a as well. Note 
that since the exact stress a is symmetric, it can be approximated to optimal accuracy by the 
symmetric subspace Eh,p- 

Next, consider the flux and trace best approximations in the quotient topology defined by ((7]) 
and (|8]). Since the exact trace it is the trace of the exact solution u, and since the exact flux a n is 
the trace of the normal components of a along each interface of df2f l , we have 

g^J^-^h^idnu) ^ \\u~n grad u\\ H i in) , 

„ inf ||<7„ - f) n ,h\\H-i/2(dn h ) ^ Ik - ^divf||if(div,fi), 

Vn,h£Qh,p 

where -Z7 gra dU G i/g(i7;V) and Tl^a £ if (div, i7;M) are suitable projections, such that their 
traces n gra< iu\E and 7TdivCrn|s on any mesh face E is in P p (E;Y) or Q P,P (E;Y). These conforming 
projectors providing approximation estimates with constants independent of p are available from 
recent works in [TTl[T2l[17l[T8 l[T9] . Specifically, QS Theorem 8.1] gives 

\\u- II siad u\\ H i( n) < C\n{p 2 ) 2 h s p^ s \u\ H s+i (n) , {s<p), (51a) 
||ct - n div cr\\ L 2( n) < C\n(p 2 ) h s p2 S \a\ H s+i^, (s<p+l). (51b) 

whenever s > 1/2 for tetrahedral meshes. The same results for cubic meshes are available from |llj . 
Note that we have chosen 77di v to be a projector into the Raviart-Thomas space. The projector 
into the Raviart-Thomas space satisfies div n dlv (j — II p div a (where 7T p denote the L 2 -orthogonal 
projection into Wh, P )- Hence, 

|| div(<7 - n dhf a)\\ L 2 (n) <Ch s P 2 S \ div a\ H s (n) , (s<p + l). (52) 

Finally, comparing the rates of convergence in (jST))) . ([51]) and ([52^1 . we find that to obtain a full 
0(h p+1 ) order of convergence, we must increase the polynomial degree of the numerical trace space 
to p + 1. Combining these observations, we have the following corollary. 

Corollary 1 (h and p convergence rates) Let fl^ be a shape regular mesh ( either tetrahedral 
or cubic) and let h denote the maximum of the diameters of its elements. Let T> and T>^ denote 
the (previously defined) discretization errors of the first and second DPG methods, resp. Using the 
spaces defined in (|49p . set 

U h = £ h , p x W h , P x M h ,p+x x F h ^ p 

for the first DPG method and 

= S h:P x W h , p x M h , p +i x F h , p x E 
for the second DPG method. Then 

v<d Hp 2 fh s p- s (\\<j\\ H ^ (n) + \\u\\ H s +1(n) ) 

£ (2) < C H \n{p 2 ) 2 h s p^{\\a\\ Hs+ , {n) + \\u\\ H s + , {n) ) 

for all 1/2 < s < p + 1. The constants Ci and Cn are independent of h and p, but dependent on 
the shape regularity and A. If Assumvtion [3T\ (isotropy) holds, then C/j is independent ofQo, so 
the second estimate does not degenerate as the Poisson ratio goes to 0.5. 
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In the same way, one can derive convergence rates for other element shapes and spaces from 
Theorems Q] and Theorem [21 

Remark 4 (Symmetric and conforming stresses) If an application demands stress approximations 
ct/j that are both symmetric and div-conforming (i.e., if one needs er^ to be in the space H (div, J?; S) 
defined in (fTJ)), then the DPG method can certainly give such approximations. We only need to 
choose a trial subspace 

£h, P c If (div, Q; S) (53) 

instead of the choice £ h , P C L 2 (f2, S) we made in above. Obviously If (div, Q; S) C L 2 (f2, M), 
so (|53|) . together with the other component spaces as set previously would result in a trial space 
Uh that is conforming in our ultraweak variational framework. Hence, Theorems [T] and [2] continue 
to apply. Notice that stability of the resulting DPG method is ensured even with this choice 
because the method adapts its test space to any given trial subspace. The first example of Eh,p 
satisfying (JSHJ) that comes to mind is the finite element of [5] . However their space is too rich because 
they had to ensure a discrete inf-sup condition. Since we have separated out the stability issue, 
we have other simpler and inexpensive options. E.g., we may choose £h,p to consist of symmetric 
matrix functions, each of whose entries are in Lh }P = {p 6 if 1 (i7,§) : p\x £ P p (K,S)} (i.e., 
continuous Lagrange finite clement functions). The locality of our test space construction is not 
destroyed with this choice. Moreover, p-optimal interpolation estimates are known for this space, 
so we can proceed as above to state an analogue of Corollary [T] Of course, the same remarks 
also apply for displacement approximations, e.g., we may choose i? 1 -conforming subspaces to 
approximate the displacement. 

7 Numerical results 

In this section, we present numerical results for the first DPG method using two test cases: a 
smooth solution on a square domain and a singular solution on an L-shaped domain. All numerical 
experiments were conducted using a pre-existing /ip-adaptive finite element package [TU] . 

7.1 Discrete spaces 

Following |13j . we considered a 2D domain Q divided into conforming or 1-irregular quadrilateral 
meshes. Let flh denote the collection of mesh elements and £ h denote the collection of mesh edges. 
A polynomial degree pk > 1 is assigned to each element and a degree pe is assigned to each mesh 
edge E. For the first method, the practical trial space is Uh = Eh x Wh x Mh x Fh where 



K h = {v : v\ K g Q PK!PK {K),VK g Q h } (54a) 

S h = (Kh) 3 (54b) 

W h = (Kh) 2 (54c) 

M h = {fi : He g P PE+ i(E),VE e £ h and fx g H^Q)} (54d) 

Fh - {ri : v\ E g P PE (E),VE g £ h } (54e) 



(and, as before, Qi^ rn is the space of bivariate polynomials which are of degree at most I in x and m 
in y). Notice that in (|54bl) . we interpret each element of (Kh) 3 as a symmetric matrix with entries 
in Kh (so the stress approximations are strongly symmetric). The edge order pe is determined 
using the maximum rule, i.e. pe is the maximum order of all elements adjacent to edge E. 

Recall that the test space is determined by the trial-to-test operator T : U i-> V, which in turn 
requires inverting the Riesz map corresponding to the inner product in the test space. In practice, 
this is solved on the discrete level by using T : U H> V which is defined as 

(fu, v) v = b(u, v) yveV (55) 



20 



J. Bramwell, L. Demkowicz, J. Gopalakrishnan, W. Qiu* 



where V is a finite dimensional subspace of V. For our implementation, V is denned as 

V={(r,v): t\ k g (Qp K ,p K ) 3 andv\ K G (Q^^.) 2 } (56) 

where = px + <5p. The default choice for the enrichment degree dp is 2. Numerical experience 
shows that this is a sufficient choice for most problems (see figure 5(d) ). 

Finally, we approximate the energy norm of the error using the error representation function 
e G V where e = T(u — Uh). Note that by the definition of T, the error representation function 
can be computed element wise by solving the variational problem 

(e, v) v = (T(u - u h ), v) v = b(u - u h , v) = l(v) - b(u h , v). (57) 

This implies that the energy norm of the error is approximated by 

\\u - u h \\ E = sup u h ,v)\ _ _ _ _ U/i )|| y = \\e\\ v . (58) 

vev \\v\\v 

From Theorem^ it is known that the energy error is equivalent to the standard norm error on U , 
so our choice of error indicator is justified assuming that the approximation of T by T is sufficient. 
In all cases, the standard test space norm is used for the inversion of the Riesz map, i.e., 

\\(r,v)\\l = \\T\\l td ^ nh) + \\v\\ 2 Hl(nh y (59) 

Note that due to the strong symmetry of functions in the discrete space Uh, the term involving 
q in (|6d[) vanishes. Hence, at the discrete level, we may (omit all g's and) work with the reduced 
test space V = H(div, fi^', §) x H (Qh;Y) instead of (|6c|) . This is why we use the norm |59|) in 
place of lf9bl). 



7.2 Test Problems 
7.2.1 Smooth Solution 

The first problem studied was a smooth solution with manufactured body force found by applying 
the elasticity equations to the exact displacements 



u x = sin(7ra;) sin(7ry) 
u y = sin(7ra;) sin(7ry) 



(60) 
(61) 



over a unit square domain S7 = (0, 1) x (0, 1) with u x ,u y prescribed on df2. For all cases, J?/j is 
initially a uniform mesh of 4 square elements. 



7.2.2 L-Shaped Steel 

The second problem considered was the classical L-shaped domain problem with the material 
properties of steel. In polar coordinates r, 9 around the re-entrant corner, the singular solution 



(62) 
(63) 
(64) 

(65) 



, [33] 


or [21 
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Smooth solution, uniform h refinement 




10" 



10 ' 



10 ' 



Smooth solution, uniform p refinement 



10' 10" 10' 

Degrees of freedom 

(a) Uniform /i-refinements. 
Fig. 1 Uniform refinement strategies for the smooth problem. 
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- Stress error 

- Displacement error 

- Energy error 



10' 10' 
Degrees of freedom 

(b) Uniform p-refinements. 



10 



where v 



2(A+m) 



is Poisson's ratio and the functions F{ff) and G(9) are given by 



F (0) = Ci sin (a + 1) 9 + C 2 cos (a + 1) 6 + C 3 sin (a - 1) 9 + C 4 cos (a - 1) i 
4 



G{9) 



a- 1 



-C 3 cos (a - 1) + C 4 sin (a - 1) 9] 



(67) 
(68) 



To determine the constants C\, C%, C^^Ca and a, we use the kinematic boundary conditions along 
the edges forming the reentrant corner (which without loss of generality we can take to be the 
edges 9 = ±37r/4). We can obtain a square integrable solution that satisfies dive = by setting 
C* 2 = C 4 = 0, C 3 = 1 and 



d = 



(l-^-)-(a+l)] sm((a-l)§7r) 



(a + l)sin((a + l)§7r) 
and letting < a < 1 be the solution of the transcendental equation 



(69) 



C\ cos 



3(a + 1)tt 



(a + 1) + cos 

+ 4 



3(a-l)vr 



1 



4 
v 



(a-1) 
'3 (a- 1)tt 



Q. 



Numerically solving for a with the material properties of steel [2T], namely A = 123 GPa, /i = 79.3 
GPa, we obtain a « 0.6038. This implies that all stress components have a singularity of strength 
(approximately) r -°- 3962 while the displacement components are smooth at the origin (but have 
singular derivatives). We can thus expect the stress components to be in (a space close to) H o eo3S ^ e 
and the displacement components in ff 1 - 6038 - 6 f or e > 0. 



7.3 Convergence rates 

The observed decrease of the error as the degrees of freedom increase is shown in Figure [T] for the 
smooth solution case and in Figure [5] for the L-shaped domain. Note that when we report the "L 2 
error" , we only consider the L 2 -norm of the errors in u and a, not the error in numerical fluxes or 
traces. 

Consider the case of the smooth solution first. If we perform uniform /i-refinements, the number 
of degrees of freedom (N) is 0(h~ 2 ). From Corollary [TJ we expect to see the error decrease by 
0(h p+1 ) for the smooth solution case, i.e., 0{N~ ijp+1 ^ 2 ) in terms of N. This is confirmed in 
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Singular solution, uniform h refinement, p = 2 Singular solution, uniform p refinement 




Degrees of freedom Degrees of freedom 



(a) Uniform h refinements. (b) Uniform p refinements. 

Fig. 2 Uniform refinement strategies for the L-shaped domain. 



Figure 1(a) Also, since both displacement and stress are infinitely smooth, they converge at the 



same rate. For uniform p-refinements, exponential convergence is observed in Figure 1(b) 



In the singular case of the L-shaped domain, Figure [2] shows the observed convergence his- 
tory for uniform refinements. Since the stress variables are in JJ - 6038 -^ we expect that the best 
approximation error for stress should decrease at rate h ' 6038 , or jv~ ' 3019 . This is in agreement 
with Figure 2(a) Additionally, since displacement is in H 1603S ~ e , one might think that its best 
approximation error should decrease more or less at rate /i 1 6038 , or Ar-o-soig However in the DPG 
method the errors for both these variables are coupled together. So, while we observe the optimal 
convergence rate for the stress variable, the convergence rate for the displacement seems to be 
somehow limited by the convergence rate of the stress. For uniform p-refinements, because we are 
considering a singular solution, the convergence rate is limited by the regularity of the solution, 



so unlike Figure 1(b) no exponential convergence is observed in Figure 2(b) 



7.4 Comparison with the weakly symmetric mixed method. 

The smooth solution problem was also implemented using the weakly symmetric mixed element 
given in 5 . This mixed method uses polynomials of degree one higher than our DPG method 
for the stress trial space. Expectedly therefore, the stress approximations given by the mixed 
method were generally observed to be superior in the L 2 norm. For the displacement however, 
both methods use the same space, so it is interesting to compare the displacement errors. This is 
done in Figure [3J The DPG method delivers lower displacement errors in the higher order case. 
In the lowest order case (not shown in the figure) the mixed method performs slightly better. 



7.5 Locking experiments 

In Figure [4] we show numerical evidence of the locking-free property of the DPG method. The 
figure shows convergence curves for various values of Poisson ratio close to the limiting value of 
0.5. We used piecewise bilinear elements with homogeneous material data. The convergence curves 
in Figure 4(a) show hardly any difference as v approaches 0.5. To be clearer, we also plot the ratio 
of the L 2 discretization error to the best approximation error (in a and u combined) in Figure [4 (b)| 
We see that the ratio remains close to the optimal value of 1.0 even as v approaches 0.5. 
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^mooth solution, uniform h refinements, order 1 displacement 



- Mixed, displacement error; 

- DPG, displacement error 




^mooth solution, uniform h refinements, order 2 displacement 

Mixed, displacement error: 
DPG, displacement error 
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Fig. 3 The DPG method vs. the mixed method. 
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Smooth solution, nearly incompressible, p = 1 
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Smooth solution, nearly incompressible, p = 1 
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(a) Convergence curves for various v 
Fig. 4 Illustration of locking-free convergence (smooth data case) 



(b) The ratio of the discretization error to the error in 
best approximation as a function of v 



7.6 Adaptivity 

All our adaptive schemes are based on the "greedy" strategy described in [13] ■ This means that 
all elements which contribute 50% of the maximum element error to the previously described 
error representation function are marked for refinement. For hp- adaptivity, we used the strategy 
suggested in pQ, i.e., if an element contains the singularity, it is /i-refined, otherwise it is p-refined. 

Figure [S] shows results from both adaptivity schemes. Note that a nearly optimal rate of 
0(N )is observed for the ft-adaptivity scheme. The /ip-adaptive scheme results in an optimal 
rate of Q(N 
urc 



1 ' 5 ). Finally, Figure 5(e) shows the hp mesh obtained after 12 iterations and Fig- 
5(f) shows one component of the corresponding solution. The group relative L 2 error is reduced 



to 0.9%. 



7.7 Approximation of optimal test functions 



Figure 5(d) shows the effect of Sp as seen in the h- adaptive process for the L-shaped domain 



problem. This measures the effect of approximating optimal test functions using the operator T 
as opposed to T, Since the curves for Sp = 2,3,4 are coincident, it appears that we are sufficiently 
approximating the optimal test functions. 
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Singular solution, adaptive h refinement, p = 2 



Singular solution, adaptive hp refinement 
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(a) Adaptive h refinements. 
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(b) Adaptive hp refinements. 
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(c) A comparison of refinement strategies. (d) Various optimal test function approximations. 
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(e) The hp mesh after 12 iterations. Element degrees 
are represented by color. (The color scale is tied to 
P + l.) 



(f) The ^-component of the computed displacement (u x ). 
The color scale indicates value of u x . 



Fig. 5 Results from the adaptive scheme for the L-shaped domain. 



A DPG method for elasticity 



25 



A A property of the weakly symmetric mixed formulation 

We consider a mixed method for linear elasticity with weakly imposed stress symmetry. The method we consider 
differs from a standard method [4] only in that it has an extra Lagrange multiplier. It is well known that the mixed 
formulation does not lock (see e.g,. l71[51l32p for homogeneous isotropic material parameters. In this appendix, we 
will provide a stability result for slightly more general materials. Note however, that the main goal of this appendix 
is to establish stability estimates for the mixed method in the form needed for the analysis of the DPG scheme in 
the earlier sections. 

The formulation reads as follows: Find (cr,u,p,a) £ ff(div, Q; M) X L 2 (i?;V) X L 2 (i?;K) X R satisfying 



(Aa, r) n + (u, divr)„ + (p, r) n + (aQg 1 AI, r) n = (F u r) n , (70a) 

(div<x,u) fi = (F 2 ,v) n , (70b) 

(°,v)n = (F 3 ,q)n, (70c) 

{a,bQ- 1 AI) n = {F 4 „bI)n, (70d) 



for all (t, v, ih b) £ H(div, Q; M), xL 2 (H; V) x L 2 (H; K) x R. Recall that Q is as defined in JT3). This formulation, 
specifically {70djl , is motivated by the same constraint that motivated the second DPG method, namely Hist . 

Theorem 5 Let (F u F 2 , F s , F 4 ) £ L 2 (S7;M) x L 2 (J?;V) X L 2 (J?;M) X L 2 (J7;M). Then: 

1. Problem (1 70 ft is uniquely solvable and the solution component u is in fact in Hg(f2;Y). 

2. There is a positive constant Co such that 

h\\ H (d,v,n) + NlHi(fi) + llPllft + M < CoOl^illft + ||F 2 || fi + ||F 3 ]| fi + ||F 4 || fi ). (71) 

3. In addition, if A ssumvt ion [371 holds, then the constant Cq in <71H takes the form 

Co = c 1 P - 1 B 4 (\\A\\ +P + 1) 2 (||A|| + B), (72) 
where c\ is a positive constant independent of A. 

We will use the Babuska-Brezzi theory [7] and results from [2] to prove this theorem. In order to verify the 
conditions of the theory, we will use the following lemma. 

Lemma 7 If Assumvtion [3l\ holds, then there is a positive constant co independent of the material coefficient A 
such that 

\\tr V \\l < cgB 2 (|| ¥ . D || 2 2 + ||div^|| 2 2 ), (73) 
for any ip £ H(div, Q; M) which satisfies 

f tr(Aip) = 0. (74) 
J n 

Proof This proof is similar to a proof in [7]. We can apply a standard regular right inverse of divergence to tr{Atp>) 
since 1741 . Hence there exists a constant co > and r\ £ Hq(/2;V) such that 

divr; = Qq 1 ti(A<p), \\v\\ Ri(fi) < c oQo 1 H tr (^v)llr2- 
By the isotropy assumption - sec (1206 - we have that tr(Aip) = Q tr ip. This implies that 

divr; = (QQq 1 )tr<^>, IMIj/i(n) < coQ^WQ^fWn- 

Then, since QQq 1 > 1 a.e., we have 

ll tri / : 'll?2 < {(QQo 1 )tri/3, tri/3) fi = (div^tr^)^ = ((div7?)J, tp)n 

= N(Wn - (Vv)D,V>)n = -N(r,,div<p) n - N((Vr]) D ,ip) a 
= -N(r), div <p)a- N(Vn, <p D )n 

<^NlHMft) (IIVDllfi + lldiv^H 2 ,) 1 / 2 

< coJVQ-^IQtr^Hn (I^dII 2 , + || div V || 2 2 .) 1/2 

Setting co = coiV, the lemma is proved. 

Proof (Proof of Theorem{^j) To apply the Babuska-Brezzi theory, we need to verify two conditions: (i) the coercivity 
on kernel, and (ii) the inf-sup condition. 

Step (i). Coercivity on kernel: Define the kernel space 

V = {t £ H(div, Q;M) : divr = 0, r' = r, / tr(Ar) = 0}. 
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Clearly, if A is uniformly coercive, then there is a positive constant ci, depending on A, such that 

c llMl!f(div,«) < (Ar,r)n, VreVb- (75) 

If in addition Assumption 1311 holds, then we can give the dependence of c\ on Po and Q, as we see now. For 
any r G V , 

(Ar,r) n > (Pr D ,T D ) n >P \\T D \\ 2 n . 
Using Lemma [7] and the fact that divT = 0, we have that 

IMIfi > (5 B)- 2 ||trr||? 2 . 

Since ||t||o = \\tjj ||^ 2 + A r_1 ||trr||^ 2 , we have that 

\\r\\l < {l + N^clB^WmWl. 

So, for any t E Vq, we have that 

-Po h ||2 _ -Po ii ||2 ^ i /. \ 

l + 7V-i S 2i?2 1|T|| H(div,fi, - i + Tv-igj^lMlfi ^ ( A ^)n 
and we conclude that (1 7511 holds with 

Cl = ™ 

in the isotropic case. 

Step (ii). Inf-sup condition: The inf-sup condition will follow once we show that there is a positive constant C2 
such that for any (u,p,a) G L 2 (fi;V) X L 2 (i?;K) X R, there is a t G H(drv,Q\M) satisfying 

(u.divr)^ + (p,r) fi + (aQg 1 AI , t)q > c 2 \\t \\ H (div,n) GM|n + \\p\\n + M) ■ (77) 

To this end, we first recall [3] Theorem 11.1]. Accordingly, there is a to G H(div, i?;M) and eg > such that 
divro = u, skwro = p, and 

IN|]H(div,0) < c 3 (||«||fi + \\p\\n)- (78) 

The constant C3 depends only on Q. 

To prove J77j , we choose r of the form r = to + \I where A G R. Obviously, 

t G H (div, n, M), divr = «, and skw r = p, 

for any A G K. So, to show the estimate 1771 . we need only choose A G K such that 

{aQ^AI^n = (aQ-'AI^o + XI) n = \a\ 2 , 

x= a-Q^(AI,r n)n 

Then, by 1781 1 . there is a positive constant C2 such that (1 77 1> holds. 

If the material is isotropic, then the dependence of C2 on the components of A can be tracked, as follows. 
Observe that since QQq > 1, we have 

Qo 1 f tr(AI) = N f Qo 1 Q>N\H\, 

\Qo 1 (AI,t )„\= f {Q- x Q)iTT <c 4 B||ro|| fi , 
J n 

with a constant C4 depending only on Q. Using this in 1791 1. we have 
MlfT(div,«) < ll T o||£r(div,f2) + 11-^ llff(div,fi) 

<, ,„ I, , || I, > , c 3 c A B{\\u\\ n + Ijpjln) + |q| 
< C3 ( tt o + P n) H . . 

Thus there is a constant C5 depending only on Q such that 

IMIfr(div,fl) < c 5 -B(|M|fi + ||p||fi + |o|), 

and therefore (1 77 6 holds with 

c 2 = (cs^r 1 (80) 

in the isotropic case. 
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Step (iii). Since we have verified the two conditions of the Babuska-Brezzi theory, we conclude that there is a 
unique solution (a,u,p,a) £ H(div, Q; M) X L 2 (0;Y) X L 2 (Q;K) X R. Moreover, the theory guarantees - see [7] 
Eq. (1.29)-(1.30)] - that the stability estimate 

lklltf(div,n) + Nlu + IIPII + kl < ^(ll-Flllfl + \\Pa\\a + \\F 3 \\(2 + \\Fi\\n) 

holds with 

C = c- 1 c- 2 (\\A\\+c 1 + c 2 f. 
If in addition, Assumption 1311 holds, then by (1761 ) and 1801 . 

C = Pq 1 {1 + c 2 B 2 )cIb 2 {\\A\\ + {c 5 B)~ 1 +P ) 2 . (81) 

Step (iv). To prove that u is in fact in H 1 (i?,V), we observe that by choosing r E Z>(i?;M) arbitrarily, we can 
conclude that the equality 

Aa - Vu + p + aQ^AI = Fi (82) 

holds in the sense of distributions. Hence u £ H 1 (i7;V). Consequently, we may integrate l !70aH by parts for any 
r e H(div, f2;M) and use (|82j to conclude that u e H^(H; V). By (gg) , 

IMIff(div,fi) + IkllHi(fi) + IIpII + M < C*(2+ ||A|| +Qo 1 ||A/|| fi )(||Fi|| fi + ||F 2 || fi + \\F 3 \\ n + ||F 4 ||r2). 

We have thus proved J7T} with Co = C(2 + ||A|| + Q~ 1 ||A7|| 

If in addition, the material is isotropic, then by H81II , the constant Co can be written as 

C = P^l + c 2 B 2 )c 2 B 2 {\\A\\ + (csB)- 1 + P ) 2 (2 + \\A\\ + N\Q\B). 

Since B > 1, we conclude that there is a positive constant ci such that Co < 5iPq — 1 B 4 (||A|| + Po + 1) 2 (||A|| + B), 
as stated in 11721 . 
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